module ObservationsCore

using Dates

"""
    ObservationRecord

Single observation record with all metadata fields.
This corresponds to a single row in an observation file.
"""
Base.@kwdef struct ObservationRecord
    # Basic observation data
    rec::String                    # receptor/station name
    yyyymmdd::Int                 # date as YYYYMMDD
    hhmmss::Int                   # time as HHMMSS
    jdate::Float64                # Julian date
    avetime::Float64              # averaging time (hours)
    conc::Float64                 # observed concentration
    err::Float64                  # measurement error
    num::Int                      # number/flag field

    # Additional metadata (can be populated later)
    datetime::Union{DateTime, Nothing} = nothing
    station_index::Union{Int, Nothing} = nothing
end

"""
    Observations

Rich observation data structure with arrays for efficient processing.
This mirrors the Fortran mod_obs.f90 obs_t type but with Julia idioms.
"""
Base.@kwdef mutable struct Observations
    # Core observation data arrays
    jdates::Vector{Float64}           # Julian dates
    stations::Vector{String}          # station/receptor names
    concentrations::Vector{Float64}   # observed concentrations
    avetime::Vector{Float64}          # averaging times
    measurement_errors::Vector{Float64} # measurement errors
    datetimes::Vector{DateTime}       # parsed datetime objects

    # Optional model comparison fields (initialized as empty, filled during inversion)
    total_errors::Vector{Float64}     # total observation errors (measurement + model)
    background::Vector{Float64}       # background concentrations
    model_values::Vector{Float64}     # modeled concentrations
    residuals::Vector{Float64}        # model - observation differences

    # Species-specific contributions (for CO2)
    nee_contrib::Vector{Float64}      # NEE contribution
    fff_contrib::Vector{Float64}      # fossil fuel contribution
    ocean_contrib::Vector{Float64}    # ocean contribution

    # Prior and posterior values
    prior_values::Vector{Float64}     # prior state contribution
    posterior_values::Vector{Float64} # posterior state contribution
    cinipos::Vector{Float64}          # initial mixing ratio contribution

    # Metadata and indices
    station_indices::Vector{Int}      # map to station/receptor list
    valid_flags::Vector{Bool}         # quality control flags

    # Lookup structures for efficient access
    by_station::Dict{String, Vector{Int}}  # indices by station name
    by_date::Dict{Date, Vector{Int}}       # indices by date

    # File metadata
    source_files::Vector{String}      # track which files data came from
    obsfiles_list::Vector{String}     # list from obsfiles.txt if used
end

"""
    empty_observations() -> Observations

Create empty Observations structure.
"""
function empty_observations()
    return Observations(
        jdates = Float64[],
        stations = String[],
        concentrations = Float64[],
        avetime = Float64[],
        measurement_errors = Float64[],
        datetimes = DateTime[],
        total_errors = Float64[],
        background = Float64[],
        model_values = Float64[],
        residuals = Float64[],
        nee_contrib = Float64[],
        fff_contrib = Float64[],
        ocean_contrib = Float64[],
        prior_values = Float64[],
        posterior_values = Float64[],
        cinipos = Float64[],
        station_indices = Int[],
        valid_flags = Bool[],
        by_station = Dict{String, Vector{Int}}(),
        by_date = Dict{Date, Vector{Int}}(),
        source_files = String[],
        obsfiles_list = String[]
    )
end

"""
    length(obs::Observations) -> Int

Return the number of observations.
"""
Base.length(obs::Observations) = length(obs.jdates)

"""
    isempty(obs::Observations) -> Bool

Check if observations structure is empty.
"""
Base.isempty(obs::Observations) = isempty(obs.jdates)

"""
    add_observation!(obs::Observations, record::ObservationRecord,
                     station_idx::Int=0, source_file::String="")

Add a single observation record to the Observations structure.
"""
function add_observation!(obs::Observations, record::ObservationRecord;
                         station_idx::Int=0, source_file::String="")
    # Add to main arrays
    push!(obs.jdates, record.jdate)
    push!(obs.stations, record.rec)
    push!(obs.concentrations, record.conc)
    push!(obs.avetime, record.avetime)
    push!(obs.measurement_errors, record.err)
    push!(obs.station_indices, station_idx)
    push!(obs.valid_flags, true)  # Default to valid

    # Parse datetime
    dt = parse_datetime(record.yyyymmdd, record.hhmmss)
    push!(obs.datetimes, dt)

    # Initialize model arrays with zeros/NaN
    push!(obs.total_errors, record.err)  # Start with measurement error
    push!(obs.background, 0.0)
    push!(obs.model_values, 0.0)
    push!(obs.residuals, 0.0)
    push!(obs.nee_contrib, 0.0)
    push!(obs.fff_contrib, 0.0)
    push!(obs.ocean_contrib, 0.0)
    push!(obs.prior_values, 0.0)
    push!(obs.posterior_values, 0.0)
    push!(obs.cinipos, 0.0)

    # Track source file
    push!(obs.source_files, source_file)

    # Update indices
    idx = length(obs.jdates)

    # Update by_station index
    if !haskey(obs.by_station, record.rec)
        obs.by_station[record.rec] = Int[]
    end
    push!(obs.by_station[record.rec], idx)

    # Update by_date index
    date_key = Date(dt)
    if !haskey(obs.by_date, date_key)
        obs.by_date[date_key] = Int[]
    end
    push!(obs.by_date[date_key], idx)

    return idx
end

"""
    parse_datetime(yyyymmdd::Int, hhmmss::Int) -> DateTime

Convert YYYYMMDD and HHMMSS integers to DateTime.
"""
function parse_datetime(yyyymmdd::Int, hhmmss::Int)
    year = yyyymmdd ÷ 10000
    month = (yyyymmdd ÷ 100) % 100
    day = yyyymmdd % 100

    hour = hhmmss ÷ 10000
    minute = (hhmmss ÷ 100) % 100
    second = hhmmss % 100

    return DateTime(year, month, day, hour, minute, second)
end

"""
    get_station_data(obs::Observations, station::String) -> Vector{Int}

Get indices of all observations for a given station.
"""
function get_station_data(obs::Observations, station::String)
    return get(obs.by_station, station, Int[])
end

"""
    get_date_data(obs::Observations, date::Date) -> Vector{Int}

Get indices of all observations for a given date.
"""
function get_date_data(obs::Observations, date::Date)
    return get(obs.by_date, date, Int[])
end

"""
    get_stations(obs::Observations) -> Vector{String}

Get list of unique station names.
"""
function get_stations(obs::Observations)
    return collect(keys(obs.by_station))
end

"""
    get_date_range(obs::Observations) -> Tuple{Date, Date}

Get the date range of observations.
"""
function get_date_range(obs::Observations)
    if isempty(obs.datetimes)
        return Date(1900, 1, 1), Date(1900, 1, 1)
    end
    dates = Date.(obs.datetimes)
    return minimum(dates), maximum(dates)
end

"""
    filter_by_date_range!(obs::Observations, start_date::Date, end_date::Date)

Filter observations to keep only those within the specified date range.
"""
function filter_by_date_range!(obs::Observations, start_date::Date, end_date::Date)
    # Find indices to keep
    keep_indices = Int[]
    for (i, dt) in enumerate(obs.datetimes)
        date = Date(dt)
        if start_date <= date <= end_date
            push!(keep_indices, i)
        end
    end

    # Filter all arrays
    obs.jdates = obs.jdates[keep_indices]
    obs.stations = obs.stations[keep_indices]
    obs.concentrations = obs.concentrations[keep_indices]
    obs.avetime = obs.avetime[keep_indices]
    obs.measurement_errors = obs.measurement_errors[keep_indices]
    obs.datetimes = obs.datetimes[keep_indices]
    obs.total_errors = obs.total_errors[keep_indices]
    obs.background = obs.background[keep_indices]
    obs.model_values = obs.model_values[keep_indices]
    obs.residuals = obs.residuals[keep_indices]
    obs.nee_contrib = obs.nee_contrib[keep_indices]
    obs.fff_contrib = obs.fff_contrib[keep_indices]
    obs.ocean_contrib = obs.ocean_contrib[keep_indices]
    obs.prior_values = obs.prior_values[keep_indices]
    obs.posterior_values = obs.posterior_values[keep_indices]
    obs.cinipos = obs.cinipos[keep_indices]
    obs.station_indices = obs.station_indices[keep_indices]
    obs.valid_flags = obs.valid_flags[keep_indices]
    obs.source_files = obs.source_files[keep_indices]

    # Rebuild indices
    rebuild_indices!(obs)
end

"""
    rebuild_indices!(obs::Observations)

Rebuild the by_station and by_date indices after filtering.
"""
function rebuild_indices!(obs::Observations)
    # Clear existing indices
    empty!(obs.by_station)
    empty!(obs.by_date)

    # Rebuild
    for (i, (station, dt)) in enumerate(zip(obs.stations, obs.datetimes))
        # by_station
        if !haskey(obs.by_station, station)
            obs.by_station[station] = Int[]
        end
        push!(obs.by_station[station], i)

        # by_date
        date_key = Date(dt)
        if !haskey(obs.by_date, date_key)
            obs.by_date[date_key] = Int[]
        end
        push!(obs.by_date[date_key], i)
    end
end

"""
    summary_stats(obs::Observations) -> NamedTuple

Get summary statistics for the observations.
"""
function summary_stats(obs::Observations)
    if isempty(obs)
        return (count=0, stations=String[], n_stations=0, date_range=(Date(1900,1,1), Date(1900,1,1)),
                concentration_range=(0.0, 0.0), mean_concentration=0.0, mean_error=0.0)
    end

    return (
        count = length(obs),
        stations = get_stations(obs),
        n_stations = length(obs.by_station),
        date_range = get_date_range(obs),
        concentration_range = (minimum(obs.concentrations), maximum(obs.concentrations)),
        mean_concentration = sum(obs.concentrations) / length(obs.concentrations),
        mean_error = sum(obs.measurement_errors) / length(obs.measurement_errors)
    )
end

"""
    apply_measurement_error_floor!(obs::Observations, min_error::Real)

Ensure all observation uncertainties respect the floor defined in the
configuration (`measerr` in the Fortran code). This mirrors the Fortran
behaviour `max(config%measerr, err)` so that Julia and Fortran cost
functions evaluate the same weights.
"""
function apply_measurement_error_floor!(obs::Observations, min_error::Real)
    min_error <= 0 && return obs

    for i in eachindex(obs.measurement_errors)
        if obs.measurement_errors[i] < min_error
            obs.measurement_errors[i] = Float64(min_error)
        end
        if obs.total_errors[i] < min_error
            obs.total_errors[i] = Float64(min_error)
        end
    end

    return obs
end

end # module
